{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cfc395fc",
   "metadata": {},
   "source": [
    "## basics"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f0539d7",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-03T12:09:36.445647Z",
     "start_time": "2023-08-03T12:09:36.428384Z"
    }
   },
   "source": [
    "- 计算 $f(x)$ w.r.t $p(x)$（distribution） 的期望\n",
    "\n",
    "$$\n",
    "\\langle f \\rangle_{p} =\\int f(x)p(x)dx\n",
    "$$\n",
    "\n",
    "- 在一些场景下，很难尤其解析解形式，所以要转向数值解，\n",
    "    - 如果可以正确地在 $p(x)$ 中进行采样，假如采样 $N$ 次，得到 $N$ 个样本：$x_1,x_2,\\cdots,x_N$\n",
    "\n",
    "$$\n",
    "\\langle f \\rangle_{p} \\approx \\frac{1}{N}\\sum_i f(x_i)_{x_i\\sim p(x)}\n",
    "$$\n",
    "\n",
    "- 现在的问题在于 $p(x)$ 很难被采样，这个时候引入 $q(x)$（也是一个概率分布，容易被采样）来间接地计算：\n",
    "    - $q(x)$：称为重要性分布（importance distribution），所谓的重要性采样，即是从这样一个简单的重要性分布中进行采样\n",
    "        - 从 $q(x)$ 采样 $N$ 次，得到 $N$ 个样本：$x_1,x_2,\\cdots,x_N$\n",
    "    - $w(x)=\\frac{p(x)}{q(x)}$：称为 importance weight\n",
    "$$\n",
    "\\begin{split}\n",
    "\\langle f \\rangle_{p}&=\\int f(x)\\frac{p(x)}{q(x)}q(x)dx\\\\\n",
    "&=\\int f(x)w(x)q(x)dx, \\qquad w(x)=\\frac{p(x)}{q(x)}\\\\\n",
    "&=\\langle f\\cdot w\\rangle_q\\\\\n",
    "&\\approx\\frac{1}{N}\\sum_i^Nf(x_i)w(x_i)_{x_i\\sim q(x)}\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4080abd",
   "metadata": {},
   "source": [
    "## 实现重要性采样"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "33f7c74d",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-03T12:46:51.628939Z",
     "start_time": "2023-08-03T12:46:51.623386Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.special import logsumexp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc8e614f",
   "metadata": {},
   "source": [
    "- 正态分布与均匀分布\n",
    "    - $f(x) = \\frac{1}{\\sigma \\sqrt{2\\pi} } e^{-\\frac{1}{2}\\left(\\frac{x-\\mu}{\\sigma}\\right)^2}$\n",
    "        - $\\log f(x)=-\\frac12\\left(\\frac{x-\\mu}{\\sigma}\\right)^2$\n",
    "    \n",
    "    - $f(x)=\\frac{1}{b-a}$\n",
    "        - $\\log f(x)=-\\log (b-a)$\n",
    "    \n",
    "- $w_i=\\frac{p(x_i)}{q(x_i)}, \\log w_i=\\log p(x_i)-\\log q(x_i)$\n",
    "    - $w_i=\\frac{\\exp(w_i)}{\\sum_j\\exp(w_j)}$\n",
    "    - $\\log w_i=w_i-\\log\\sum_j\\exp(w_j)$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8be2d476",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-03T12:32:25.095379Z",
     "start_time": "2023-08-03T12:32:25.083243Z"
    }
   },
   "outputs": [],
   "source": [
    "class Pdf:\n",
    "    def __call__(self, x):\n",
    "        '''log prob'''\n",
    "        pass\n",
    "    def sample(self, n):\n",
    "        pass\n",
    "\n",
    "# 待估计的分布\n",
    "class Norm(Pdf):\n",
    "    def __init__(self, mu=0, sigma=1):\n",
    "        self.mu = mu\n",
    "        self.sigma = sigma\n",
    "        \n",
    "    def __call__(self, x):\n",
    "        return -0.5*(x-self.mu)**2/self.sigma**2\n",
    "    \n",
    "    # 不会被调用\n",
    "    def sample(self, n):\n",
    "        return np.random.normal(self.mu, self.sigma, n)\n",
    "\n",
    "# 重要性分布\n",
    "class Uniform(Pdf):\n",
    "    def __init__(self, low, high):\n",
    "        self.low = low\n",
    "        self.high = high\n",
    "    def __call__(self, x):\n",
    "        return np.repeat(-np.log(self.high - self.low), len(x))\n",
    "    def sample(self, n):\n",
    "        return np.random.uniform(self.low, self.high, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "af07a9ca",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-03T12:47:17.550357Z",
     "start_time": "2023-08-03T12:47:17.540590Z"
    }
   },
   "outputs": [],
   "source": [
    "class ImportanceSampler:\n",
    "    def __init__(self, p_dist, q_dist):\n",
    "        self.p_dist = p_dist\n",
    "        self.q_dist = q_dist\n",
    "        \n",
    "    def sample(self, n):\n",
    "        samples = self.q_dist.sample(n)\n",
    "        weights = self.calc_weights(samples)\n",
    "        norm_weights = weights - logsumexp(weights)\n",
    "        return samples, norm_weights\n",
    "    def calc_weights(self, samples):\n",
    "        # log(p(x)/q(x)) = log(p) - log(q)\n",
    "        return self.p_dist(samples) - self.q_dist(samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "725f05cd",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-03T13:07:17.003554Z",
     "start_time": "2023-08-03T13:07:16.993421Z"
    }
   },
   "outputs": [],
   "source": [
    "N = 200000\n",
    "target_p = Norm()\n",
    "imp_q = Uniform(-20, 30)\n",
    "sampler = ImportanceSampler(target_p, imp_q)\n",
    "biased_samples, logws = sampler.sample(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "8ce915e3",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-03T13:07:18.274353Z",
     "start_time": "2023-08-03T13:07:18.167671Z"
    }
   },
   "outputs": [],
   "source": [
    "samples = np.random.choice(biased_samples, N, p=np.exp(logws))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "882821a9",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-03T13:07:20.093651Z",
     "start_time": "2023-08-03T13:07:19.835029Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAGdCAYAAAAbudkLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAvHElEQVR4nO3df3BV9Z3/8dc1kAtkk1NCvLm5Y0DWwSwYVAg2BNwKggkpIYu6BRvnrmxpqCOS5ptkVOxsi51KrIB0RgaKDAWFaJwugu6GZgl1BRkIYGiqAWTRggRJCEJyL2HpTYzn+0e/nK+XIBJMvOTD8zFzZnI+533PfZ876n35OT+uy7ZtWwAAAAa6IdINAAAA9BSCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWH0i3UAkffHFFzpx4oRiY2Plcrki3Q4AALgCtm3r7Nmz8vl8uuGGy8/ZXNdB58SJE0pOTo50GwAA4CrU19frpptuumzNdR10YmNjJf3tg4qLi4twNwAA4EoEg0ElJyc73+OXc10HnQunq+Li4gg6AAD0Mldy2QkXIwMAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsboUdEpLS3XXXXcpNjZWHo9H06dP16FDh8JqbNvWggUL5PP51L9/f02YMEH79+8PqwmFQpo3b54SEhIUExOj3NxcHT9+PKymublZfr9flmXJsiz5/X61tLSE1Rw7dkzTpk1TTEyMEhISVFBQoLa2tq4cEgAAMFiXgs62bds0d+5cVVdXq6qqSp9//rkyMzN17tw5p+b555/XCy+8oGXLlmnv3r3yer267777dPbsWaemsLBQGzduVHl5uXbs2KHW1lbl5OSoo6PDqcnLy1Ntba0qKytVWVmp2tpa+f1+Z3tHR4emTp2qc+fOaceOHSovL9eGDRtUXFz8TT4PAABgEvsbaGpqsiXZ27Zts23btr/44gvb6/Xazz33nFPz17/+1bYsy/7tb39r27Ztt7S02H379rXLy8udmk8//dS+4YYb7MrKStu2bfvAgQO2JLu6utqp2bVrly3J/vDDD23btu3NmzfbN9xwg/3pp586Na+99prtdrvtQCBwRf0HAgFb0hXXAwCAyOvK9/c3ukYnEAhIkuLj4yVJR44cUWNjozIzM50at9ute+65Rzt37pQk1dTUqL29PazG5/MpNTXVqdm1a5csy1J6erpTM3bsWFmWFVaTmpoqn8/n1GRlZSkUCqmmpuaS/YZCIQWDwbAFAACYq8/VvtC2bRUVFenuu+9WamqqJKmxsVGSlJiYGFabmJioTz75xKmJjo7WwIEDO9VceH1jY6M8Hk+n9/R4PGE1F7/PwIEDFR0d7dRcrLS0VM8880xXDxXANeTmpyquuPboc1N7sBMAvcFVz+g8/vjjev/99/Xaa6912uZyucLWbdvuNHaxi2suVX81NV82f/58BQIBZ6mvr79sTwAAoHe7qqAzb948vfXWW/rv//5v3XTTTc641+uVpE4zKk1NTc7si9frVVtbm5qbmy9bc/LkyU7ve+rUqbCai9+nublZ7e3tnWZ6LnC73YqLiwtbAACAubp06sq2bc2bN08bN27UO++8o6FDh4ZtHzp0qLxer6qqqjRq1ChJUltbm7Zt26Zf//rXkqS0tDT17dtXVVVVmjFjhiSpoaFBdXV1ev755yVJGRkZCgQC2rNnj7773e9Kknbv3q1AIKBx48Y5Nc8++6waGhqUlJQkSdqyZYvcbrfS0tKu9vMAYBBOcwHoUtCZO3euXn31Vb355puKjY11ZlQsy1L//v3lcrlUWFiohQsXatiwYRo2bJgWLlyoAQMGKC8vz6mdPXu2iouLNWjQIMXHx6ukpEQjR47U5MmTJUnDhw/XlClTlJ+fr5UrV0qS5syZo5ycHKWkpEiSMjMzNWLECPn9fi1atEhnzpxRSUmJ8vPzmakBAACSuhh0VqxYIUmaMGFC2PiaNWs0a9YsSdITTzyh8+fP67HHHlNzc7PS09O1ZcsWxcbGOvVLly5Vnz59NGPGDJ0/f16TJk3S2rVrFRUV5dSUlZWpoKDAuTsrNzdXy5Ytc7ZHRUWpoqJCjz32mMaPH6/+/fsrLy9Pixcv7tIHAAAAzOWybduOdBOREgwGZVmWAoEAs0BAL9GV01FdwakroPfoyvc3v3UFAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADAWQQcAABirS09GBgBTdfVBhDxgEOgdmNEBAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADAWQQcAABiLoAMAAIxF0AEAAMbiycgAcBW68iRlnqIMRA4zOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADG6nLQ2b59u6ZNmyafzyeXy6VNmzaFbXe5XJdcFi1a5NRMmDCh0/aHHnoobD/Nzc3y+/2yLEuWZcnv96ulpSWs5tixY5o2bZpiYmKUkJCggoICtbW1dfWQAACAobocdM6dO6c77rhDy5Ytu+T2hoaGsOV3v/udXC6XHnzwwbC6/Pz8sLqVK1eGbc/Ly1Ntba0qKytVWVmp2tpa+f1+Z3tHR4emTp2qc+fOaceOHSovL9eGDRtUXFzc1UMCAACG6tPVF2RnZys7O/srt3u93rD1N998UxMnTtTf//3fh40PGDCgU+0FBw8eVGVlpaqrq5Weni5JWrVqlTIyMnTo0CGlpKRoy5YtOnDggOrr6+Xz+SRJS5Ys0axZs/Tss88qLi6uq4cGAAAM06PX6Jw8eVIVFRWaPXt2p21lZWVKSEjQbbfdppKSEp09e9bZtmvXLlmW5YQcSRo7dqwsy9LOnTudmtTUVCfkSFJWVpZCoZBqamou2U8oFFIwGAxbAACAubo8o9MVL7/8smJjY/XAAw+EjT/88MMaOnSovF6v6urqNH/+fP35z39WVVWVJKmxsVEej6fT/jwejxobG52axMTEsO0DBw5UdHS0U3Ox0tJSPfPMM91xaAAAoBfo0aDzu9/9Tg8//LD69esXNp6fn+/8nZqaqmHDhmnMmDHat2+fRo8eLelvFzVfzLbtsPErqfmy+fPnq6ioyFkPBoNKTk7u2kEBAIBeo8dOXb377rs6dOiQfvzjH39t7ejRo9W3b18dPnxY0t+u8zl58mSnulOnTjmzOF6vt9PMTXNzs9rb2zvN9FzgdrsVFxcXtgAAAHP1WNBZvXq10tLSdMcdd3xt7f79+9Xe3q6kpCRJUkZGhgKBgPbs2ePU7N69W4FAQOPGjXNq6urq1NDQ4NRs2bJFbrdbaWlp3Xw0AACgN+ryqavW1lZ99NFHzvqRI0dUW1ur+Ph4DR48WNLfTgn9/ve/15IlSzq9/uOPP1ZZWZm+//3vKyEhQQcOHFBxcbFGjRql8ePHS5KGDx+uKVOmKD8/37ntfM6cOcrJyVFKSookKTMzUyNGjJDf79eiRYt05swZlZSUKD8/n5kaAAAg6SpmdN577z2NGjVKo0aNkiQVFRVp1KhR+vnPf+7UlJeXy7Zt/fCHP+z0+ujoaP3xj39UVlaWUlJSVFBQoMzMTG3dulVRUVFOXVlZmUaOHKnMzExlZmbq9ttv17p165ztUVFRqqioUL9+/TR+/HjNmDFD06dP1+LFi7t6SAAAwFAu27btSDcRKcFgUJZlKRAIMAsE9BI3P1UR6Ra67OhzUyPdAmCUrnx/81tXAADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABj9eiPegLAleiNz8YB0DswowMAAIxF0AEAAMYi6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADAWQQcAABiLoAMAAIzV5aCzfft2TZs2TT6fTy6XS5s2bQrbPmvWLLlcrrBl7NixYTWhUEjz5s1TQkKCYmJilJubq+PHj4fVNDc3y+/3y7IsWZYlv9+vlpaWsJpjx45p2rRpiomJUUJCggoKCtTW1tbVQwIAAIbqctA5d+6c7rjjDi1btuwra6ZMmaKGhgZn2bx5c9j2wsJCbdy4UeXl5dqxY4daW1uVk5Ojjo4OpyYvL0+1tbWqrKxUZWWlamtr5ff7ne0dHR2aOnWqzp07px07dqi8vFwbNmxQcXFxVw8JAAAYqk9XX5Cdna3s7OzL1rjdbnm93ktuCwQCWr16tdatW6fJkydLktavX6/k5GRt3bpVWVlZOnjwoCorK1VdXa309HRJ0qpVq5SRkaFDhw4pJSVFW7Zs0YEDB1RfXy+fzydJWrJkiWbNmqVnn31WcXFxXT00AABgmB65Ruedd96Rx+PRrbfeqvz8fDU1NTnbampq1N7erszMTGfM5/MpNTVVO3fulCTt2rVLlmU5IUeSxo4dK8uywmpSU1OdkCNJWVlZCoVCqqmpuWRfoVBIwWAwbAEAAObq9qCTnZ2tsrIyvf3221qyZIn27t2re++9V6FQSJLU2Nio6OhoDRw4MOx1iYmJamxsdGo8Hk+nfXs8nrCaxMTEsO0DBw5UdHS0U3Ox0tJS55ofy7KUnJz8jY8XAABcu7p86urrzJw50/k7NTVVY8aM0ZAhQ1RRUaEHHnjgK19n27ZcLpez/uW/v0nNl82fP19FRUXOejAYJOwAAGCwHr+9PCkpSUOGDNHhw4clSV6vV21tbWpubg6ra2pqcmZovF6vTp482Wlfp06dCqu5eOamublZ7e3tnWZ6LnC73YqLiwtbAACAuXo86Jw+fVr19fVKSkqSJKWlpalv376qqqpyahoaGlRXV6dx48ZJkjIyMhQIBLRnzx6nZvfu3QoEAmE1dXV1amhocGq2bNkit9uttLS0nj4sAADQC3T51FVra6s++ugjZ/3IkSOqra1VfHy84uPjtWDBAj344INKSkrS0aNH9fTTTyshIUH333+/JMmyLM2ePVvFxcUaNGiQ4uPjVVJSopEjRzp3YQ0fPlxTpkxRfn6+Vq5cKUmaM2eOcnJylJKSIknKzMzUiBEj5Pf7tWjRIp05c0YlJSXKz89npgbANeXmpyquuPboc1N7sBPg+tPloPPee+9p4sSJzvqFa14eeeQRrVixQh988IFeeeUVtbS0KCkpSRMnTtTrr7+u2NhY5zVLly5Vnz59NGPGDJ0/f16TJk3S2rVrFRUV5dSUlZWpoKDAuTsrNzc37Nk9UVFRqqio0GOPPabx48erf//+ysvL0+LFi7v+KQAAACO5bNu2I91EpASDQVmWpUAgwCwQEEFdmfEwHTM6wNfryvc3v3UFAACMRdABAADGIugAAABjEXQAAICxuv3JyADAxcUArhXM6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADAWQQcAABiLoAMAAIxF0AEAAMYi6AAAAGP1iXQDAID/7+anKq649uhzU3uwE8AMzOgAAABjEXQAAICxuhx0tm/frmnTpsnn88nlcmnTpk3Otvb2dj355JMaOXKkYmJi5PP59C//8i86ceJE2D4mTJggl8sVtjz00ENhNc3NzfL7/bIsS5Zlye/3q6WlJazm2LFjmjZtmmJiYpSQkKCCggK1tbV19ZAAAIChuhx0zp07pzvuuEPLli3rtO1///d/tW/fPv3bv/2b9u3bpzfeeEP/8z//o9zc3E61+fn5amhocJaVK1eGbc/Ly1Ntba0qKytVWVmp2tpa+f1+Z3tHR4emTp2qc+fOaceOHSovL9eGDRtUXFzc1UMCAACG6vLFyNnZ2crOzr7kNsuyVFVVFTb24osv6rvf/a6OHTumwYMHO+MDBgyQ1+u95H4OHjyoyspKVVdXKz09XZK0atUqZWRk6NChQ0pJSdGWLVt04MAB1dfXy+fzSZKWLFmiWbNm6dlnn1VcXFxXDw0AABimx6/RCQQCcrlc+s53vhM2XlZWpoSEBN12220qKSnR2bNnnW27du2SZVlOyJGksWPHyrIs7dy506lJTU11Qo4kZWVlKRQKqaam5pK9hEIhBYPBsAUAAJirR28v/+tf/6qnnnpKeXl5YTMsDz/8sIYOHSqv16u6ujrNnz9ff/7zn53ZoMbGRnk8nk7783g8amxsdGoSExPDtg8cOFDR0dFOzcVKS0v1zDPPdNfhAQCAa1yPBZ329nY99NBD+uKLL7R8+fKwbfn5+c7fqampGjZsmMaMGaN9+/Zp9OjRkiSXy9Vpn7Zth41fSc2XzZ8/X0VFRc56MBhUcnJy1w4MAAD0Gj1y6qq9vV0zZszQkSNHVFVV9bXXy4wePVp9+/bV4cOHJUler1cnT57sVHfq1ClnFsfr9XaauWlublZ7e3unmZ4L3G634uLiwhYAAGCubg86F0LO4cOHtXXrVg0aNOhrX7N//361t7crKSlJkpSRkaFAIKA9e/Y4Nbt371YgENC4ceOcmrq6OjU0NDg1W7ZskdvtVlpaWjcfFQAA6I26fOqqtbVVH330kbN+5MgR1dbWKj4+Xj6fT//8z/+sffv26T//8z/V0dHhzLrEx8crOjpaH3/8scrKyvT9739fCQkJOnDggIqLizVq1CiNHz9ekjR8+HBNmTJF+fn5zm3nc+bMUU5OjlJSUiRJmZmZGjFihPx+vxYtWqQzZ86opKRE+fn5zNQAAABJVzGj895772nUqFEaNWqUJKmoqEijRo3Sz3/+cx0/flxvvfWWjh8/rjvvvFNJSUnOcuFuqejoaP3xj39UVlaWUlJSVFBQoMzMTG3dulVRUVHO+5SVlWnkyJHKzMxUZmambr/9dq1bt87ZHhUVpYqKCvXr10/jx4/XjBkzNH36dC1evPibfiYAAMAQLtu27Ug3ESnBYFCWZSkQCDALBHSjrvwwJa4eP+qJ61VXvr/5rSsAAGAsgg4AADAWQQcAABiLoAMAAIxF0AEAAMYi6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABj9Yl0AwB6h5ufqoh0CwDQZczoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxuhx0tm/frmnTpsnn88nlcmnTpk1h223b1oIFC+Tz+dS/f39NmDBB+/fvD6sJhUKaN2+eEhISFBMTo9zcXB0/fjysprm5WX6/X5ZlybIs+f1+tbS0hNUcO3ZM06ZNU0xMjBISElRQUKC2trauHhIA9Eo3P1XRpQW4HnU56Jw7d0533HGHli1bdsntzz//vF544QUtW7ZMe/fuldfr1X333aezZ886NYWFhdq4caPKy8u1Y8cOtba2KicnRx0dHU5NXl6eamtrVVlZqcrKStXW1srv9zvbOzo6NHXqVJ07d047duxQeXm5NmzYoOLi4q4eEgAAMJTLtm37ql/scmnjxo2aPn26pL/N5vh8PhUWFurJJ5+U9LfZm8TERP3617/WT37yEwUCAd14441at26dZs6cKUk6ceKEkpOTtXnzZmVlZengwYMaMWKEqqurlZ6eLkmqrq5WRkaGPvzwQ6WkpOgPf/iDcnJyVF9fL5/PJ0kqLy/XrFmz1NTUpLi4uK/tPxgMyrIsBQKBK6oHrmfMCPR+R5+bGukWgG7Rle/vbr1G58iRI2psbFRmZqYz5na7dc8992jnzp2SpJqaGrW3t4fV+Hw+paamOjW7du2SZVlOyJGksWPHyrKssJrU1FQn5EhSVlaWQqGQampqLtlfKBRSMBgMWwAAgLm6Neg0NjZKkhITE8PGExMTnW2NjY2Kjo7WwIEDL1vj8Xg67d/j8YTVXPw+AwcOVHR0tFNzsdLSUueaH8uylJycfBVHCQAAeoseuevK5XKFrdu23WnsYhfXXKr+amq+bP78+QoEAs5SX19/2Z4AAEDv1q1Bx+v1SlKnGZWmpiZn9sXr9aqtrU3Nzc2XrTl58mSn/Z86dSqs5uL3aW5uVnt7e6eZngvcbrfi4uLCFgAAYK5uDTpDhw6V1+tVVVWVM9bW1qZt27Zp3LhxkqS0tDT17ds3rKahoUF1dXVOTUZGhgKBgPbs2ePU7N69W4FAIKymrq5ODQ0NTs2WLVvkdruVlpbWnYcFAAB6qT5dfUFra6s++ugjZ/3IkSOqra1VfHy8Bg8erMLCQi1cuFDDhg3TsGHDtHDhQg0YMEB5eXmSJMuyNHv2bBUXF2vQoEGKj49XSUmJRo4cqcmTJ0uShg8frilTpig/P18rV66UJM2ZM0c5OTlKSUmRJGVmZmrEiBHy+/1atGiRzpw5o5KSEuXn5zNTAwAAJF1F0Hnvvfc0ceJEZ72oqEiS9Mgjj2jt2rV64okndP78eT322GNqbm5Wenq6tmzZotjYWOc1S5cuVZ8+fTRjxgydP39ekyZN0tq1axUVFeXUlJWVqaCgwLk7Kzc3N+zZPVFRUaqoqNBjjz2m8ePHq3///srLy9PixYu7/ikAAAAjfaPn6PR2PEcHuHI8R6f34zk6MEXEnqMDAABwLSHoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADAWQQcAABiLoAMAAIxF0AEAAMYi6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsfpEugEAkXPzUxWRbgEAehQzOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjNXtQefmm2+Wy+XqtMydO1eSNGvWrE7bxo4dG7aPUCikefPmKSEhQTExMcrNzdXx48fDapqbm+X3+2VZlizLkt/vV0tLS3cfDgAA6MW6Pejs3btXDQ0NzlJVVSVJ+sEPfuDUTJkyJaxm8+bNYfsoLCzUxo0bVV5erh07dqi1tVU5OTnq6OhwavLy8lRbW6vKykpVVlaqtrZWfr+/uw8HAAD0Yt3+HJ0bb7wxbP25557TLbfconvuuccZc7vd8nq9l3x9IBDQ6tWrtW7dOk2ePFmStH79eiUnJ2vr1q3KysrSwYMHVVlZqerqaqWnp0uSVq1apYyMDB06dEgpKSndfVgAAKAX6tFrdNra2rR+/Xr96Ec/ksvlcsbfeecdeTwe3XrrrcrPz1dTU5OzraamRu3t7crMzHTGfD6fUlNTtXPnTknSrl27ZFmWE3IkaezYsbIsy6m5lFAopGAwGLYAAABz9WjQ2bRpk1paWjRr1ixnLDs7W2VlZXr77be1ZMkS7d27V/fee69CoZAkqbGxUdHR0Ro4cGDYvhITE9XY2OjUeDyeTu/n8XicmkspLS11rumxLEvJycndcJQAAOBa1aM/AbF69WplZ2fL5/M5YzNnznT+Tk1N1ZgxYzRkyBBVVFTogQce+Mp92bYdNiv05b+/quZi8+fPV1FRkbMeDAYJOwAAGKzHgs4nn3yirVu36o033rhsXVJSkoYMGaLDhw9Lkrxer9ra2tTc3Bw2q9PU1KRx48Y5NSdPnuy0r1OnTikxMfEr38vtdsvtdl/N4QAAgF6ox05drVmzRh6PR1OnTr1s3enTp1VfX6+kpCRJUlpamvr27evcrSVJDQ0Nqqurc4JORkaGAoGA9uzZ49Ts3r1bgUDAqQEAAOiRGZ0vvvhCa9as0SOPPKI+ff7/W7S2tmrBggV68MEHlZSUpKNHj+rpp59WQkKC7r//fkmSZVmaPXu2iouLNWjQIMXHx6ukpEQjR4507sIaPny4pkyZovz8fK1cuVKSNGfOHOXk5HDHFQAAcPRI0Nm6dauOHTumH/3oR2HjUVFR+uCDD/TKK6+opaVFSUlJmjhxol5//XXFxsY6dUuXLlWfPn00Y8YMnT9/XpMmTdLatWsVFRXl1JSVlamgoMC5Oys3N1fLli3ricMBAAC9lMu2bTvSTURKMBiUZVkKBAKKi4uLdDvAt+7mpyoi3QK+RUefu/ylBEBv0ZXvb37rCgAAGIugAwAAjEXQAQAAxurRBwYCAK4dXbkmi+t5YApmdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsXgyMmAQfo0cAMIxowMAAIxF0AEAAMYi6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjdXvQWbBggVwuV9ji9Xqd7bZta8GCBfL5fOrfv78mTJig/fv3h+0jFApp3rx5SkhIUExMjHJzc3X8+PGwmubmZvn9flmWJcuy5Pf71dLS0t2HAwAAerEemdG57bbb1NDQ4CwffPCBs+3555/XCy+8oGXLlmnv3r3yer267777dPbsWaemsLBQGzduVHl5uXbs2KHW1lbl5OSoo6PDqcnLy1Ntba0qKytVWVmp2tpa+f3+njgcAADQS/XpkZ326RM2i3OBbdv6zW9+o5/97Gd64IEHJEkvv/yyEhMT9eqrr+onP/mJAoGAVq9erXXr1mny5MmSpPXr1ys5OVlbt25VVlaWDh48qMrKSlVXVys9PV2StGrVKmVkZOjQoUNKSUnpicMCAAC9TI/M6Bw+fFg+n09Dhw7VQw89pL/85S+SpCNHjqixsVGZmZlOrdvt1j333KOdO3dKkmpqatTe3h5W4/P5lJqa6tTs2rVLlmU5IUeSxo4dK8uynJpLCYVCCgaDYQsAADBXtwed9PR0vfLKK/qv//ovrVq1So2NjRo3bpxOnz6txsZGSVJiYmLYaxITE51tjY2Nio6O1sCBAy9b4/F4Or23x+Nxai6ltLTUuabHsiwlJyd/o2MFAADXtm4POtnZ2XrwwQc1cuRITZ48WRUVFZL+dorqApfLFfYa27Y7jV3s4ppL1X/dfubPn69AIOAs9fX1V3RMAACgd+qRa3S+LCYmRiNHjtThw4c1ffp0SX+bkUlKSnJqmpqanFker9ertrY2NTc3h83qNDU1ady4cU7NyZMnO73XqVOnOs0WfZnb7Zbb7e6OwwIAo938VMUV1x59bmoPdgJ8Mz3+HJ1QKKSDBw8qKSlJQ4cOldfrVVVVlbO9ra1N27Ztc0JMWlqa+vbtG1bT0NCguro6pyYjI0OBQEB79uxxanbv3q1AIODUAAAAdPuMTklJiaZNm6bBgwerqalJv/rVrxQMBvXII4/I5XKpsLBQCxcu1LBhwzRs2DAtXLhQAwYMUF5eniTJsizNnj1bxcXFGjRokOLj41VSUuKcCpOk4cOHa8qUKcrPz9fKlSslSXPmzFFOTg53XAEAAEe3B53jx4/rhz/8oT777DPdeOONGjt2rKqrqzVkyBBJ0hNPPKHz58/rscceU3Nzs9LT07VlyxbFxsY6+1i6dKn69OmjGTNm6Pz585o0aZLWrl2rqKgop6asrEwFBQXO3Vm5ublatmxZdx8OAADoxVy2bduRbiJSgsGgLMtSIBBQXFxcpNsBvrGuXFcBdBeu0cG3rSvf3/zWFQAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICx+kS6AQCXd/NTFZFuAQB6LWZ0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADAWQQcAABiLXy8HAHwjNz9V0aX6o89N7aFOgM6Y0QEAAMYi6AAAAGMRdAAAgLEIOgAAwFjdHnRKS0t11113KTY2Vh6PR9OnT9ehQ4fCambNmiWXyxW2jB07NqwmFApp3rx5SkhIUExMjHJzc3X8+PGwmubmZvn9flmWJcuy5Pf71dLS0t2HBAAAeqluDzrbtm3T3LlzVV1draqqKn3++efKzMzUuXPnwuqmTJmihoYGZ9m8eXPY9sLCQm3cuFHl5eXasWOHWltblZOTo46ODqcmLy9PtbW1qqysVGVlpWpra+X3+7v7kAAAQC/V7beXV1ZWhq2vWbNGHo9HNTU1+t73vueMu91ueb3eS+4jEAho9erVWrdunSZPnixJWr9+vZKTk7V161ZlZWXp4MGDqqysVHV1tdLT0yVJq1atUkZGhg4dOqSUlJTuPjQAANDL9Pg1OoFAQJIUHx8fNv7OO+/I4/Ho1ltvVX5+vpqampxtNTU1am9vV2ZmpjPm8/mUmpqqnTt3SpJ27doly7KckCNJY8eOlWVZTs3FQqGQgsFg2AIAAMzVo0HHtm0VFRXp7rvvVmpqqjOenZ2tsrIyvf3221qyZIn27t2re++9V6FQSJLU2Nio6OhoDRw4MGx/iYmJamxsdGo8Hk+n9/R4PE7NxUpLS53reSzLUnJycncdKgAAuAb16JORH3/8cb3//vvasWNH2PjMmTOdv1NTUzVmzBgNGTJEFRUVeuCBB75yf7Zty+VyOetf/vurar5s/vz5KioqctaDwSBhBwAAg/VY0Jk3b57eeustbd++XTfddNNla5OSkjRkyBAdPnxYkuT1etXW1qbm5uawWZ2mpiaNGzfOqTl58mSnfZ06dUqJiYmXfB+32y232321hwR0m64+Mh8AcHW6/dSVbdt6/PHH9cYbb+jtt9/W0KFDv/Y1p0+fVn19vZKSkiRJaWlp6tu3r6qqqpyahoYG1dXVOUEnIyNDgUBAe/bscWp2796tQCDg1AAAgOtbt8/ozJ07V6+++qrefPNNxcbGOtfLWJal/v37q7W1VQsWLNCDDz6opKQkHT16VE8//bQSEhJ0//33O7WzZ89WcXGxBg0apPj4eJWUlGjkyJHOXVjDhw/XlClTlJ+fr5UrV0qS5syZo5ycHO64AgAAknog6KxYsUKSNGHChLDxNWvWaNasWYqKitIHH3ygV155RS0tLUpKStLEiRP1+uuvKzY21qlfunSp+vTpoxkzZuj8+fOaNGmS1q5dq6ioKKemrKxMBQUFzt1Zubm5WrZsWXcfEgAA6KVctm3bkW4iUoLBoCzLUiAQUFxcXKTbwXWEa3RwPTv63NRIt4Berivf3/zWFQAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsXr0t64AALhYVx6vwK3o+KaY0QEAAMYi6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIvn6ADdpCvPBgEAfDuY0QEAAMYi6AAAAGNx6goAcM3i5yLwTTGjAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFs/RAb4CP+kAAL0fMzoAAMBYzOgAAIzAU5RxKczoAAAAYxF0AACAsQg6AADAWFyjg+sKd1IBwPWFoAMAuO509X96uHi59+r1p66WL1+uoUOHql+/fkpLS9O7774b6ZYAAMA1olfP6Lz++usqLCzU8uXLNX78eK1cuVLZ2dk6cOCABg8eHOn28C3hdBQA4Ku4bNu2I93E1UpPT9fo0aO1YsUKZ2z48OGaPn26SktLv/b1wWBQlmUpEAgoLi6uJ1tFDyLoALiWcJqr53Xl+7vXzui0tbWppqZGTz31VNh4Zmamdu7cecnXhEIhhUIhZz0QCEj62weGa0vqL/4r0i0AwFUZ/H9+f8W1dc9k9WAn5rrwvX0lczW9Nuh89tln6ujoUGJiYth4YmKiGhsbL/ma0tJSPfPMM53Gk5OTe6RHAAAux/pNpDvo3c6ePSvLsi5b02uDzgUulyts3bbtTmMXzJ8/X0VFRc76F198oTNnzmjQoEFf+RqECwaDSk5OVn19Paf7ehCf87eDz/nbwef87biePmfbtnX27Fn5fL6vre21QSchIUFRUVGdZm+ampo6zfJc4Ha75Xa7w8a+853v9FSLRouLizP+X6RrAZ/zt4PP+dvB5/ztuF4+56+bybmg195eHh0drbS0NFVVVYWNV1VVady4cRHqCgAAXEt67YyOJBUVFcnv92vMmDHKyMjQSy+9pGPHjunRRx+NdGsAAOAa0KuDzsyZM3X69Gn98pe/VENDg1JTU7V582YNGTIk0q0Zy+126xe/+EWnU4DoXnzO3w4+528Hn/O3g8/50nr1c3QAAAAup9deowMAAPB1CDoAAMBYBB0AAGAsgg4AADAWQQffWCgU0p133imXy6Xa2tpIt2OUo0ePavbs2Ro6dKj69++vW265Rb/4xS/U1tYW6dZ6veXLl2vo0KHq16+f0tLS9O6770a6JaOUlpbqrrvuUmxsrDwej6ZPn65Dhw5Fui3jlZaWyuVyqbCwMNKtXDMIOvjGnnjiiSt6DDe67sMPP9QXX3yhlStXav/+/Vq6dKl++9vf6umnn450a73a66+/rsLCQv3sZz/Tn/70J/3jP/6jsrOzdezYsUi3Zoxt27Zp7ty5qq6uVlVVlT7//HNlZmbq3LlzkW7NWHv37tVLL72k22+/PdKtXFO4vRzfyB/+8AcVFRVpw4YNuu222/SnP/1Jd955Z6TbMtqiRYu0YsUK/eUvf4l0K71Wenq6Ro8erRUrVjhjw4cP1/Tp01VaWhrBzsx16tQpeTwebdu2Td/73vci3Y5xWltbNXr0aC1fvly/+tWvdOedd+o3v/lNpNu6JjCjg6t28uRJ5efna926dRowYECk27luBAIBxcfHR7qNXqutrU01NTXKzMwMG8/MzNTOnTsj1JX5AoGAJPHPbg+ZO3eupk6dqsmTJ0e6lWtOr34yMiLHtm3NmjVLjz76qMaMGaOjR49GuqXrwscff6wXX3xRS5YsiXQrvdZnn32mjo6OTj/+m5iY2OlHgtE9bNtWUVGR7r77bqWmpka6HeOUl5dr37592rt3b6RbuSYxo4MwCxYskMvluuzy3nvv6cUXX1QwGNT8+fMj3XKvdKWf85edOHFCU6ZM0Q9+8AP9+Mc/jlDn5nC5XGHrtm13GkP3ePzxx/X+++/rtddei3Qrxqmvr9dPf/pTrV+/Xv369Yt0O9ckrtFBmM8++0yfffbZZWtuvvlmPfTQQ/qP//iPsC+Gjo4ORUVF6eGHH9bLL7/c0632alf6OV/4D9eJEyc0ceJEpaena+3atbrhBv4f5Wq1tbVpwIAB+v3vf6/777/fGf/pT3+q2tpabdu2LYLdmWfevHnatGmTtm/frqFDh0a6HeNs2rRJ999/v6Kiopyxjo4OuVwu3XDDDQqFQmHbrkcEHVyVY8eOKRgMOusnTpxQVlaW/v3f/13p6em66aabItidWT799FNNnDhRaWlpWr9+/XX/H63ukJ6errS0NC1fvtwZGzFihP7pn/6Ji5G7iW3bmjdvnjZu3Kh33nlHw4YNi3RLRjp79qw++eSTsLF//dd/1T/8wz/oySef5FShuEYHV2nw4MFh63/3d38nSbrlllsIOd3oxIkTmjBhggYPHqzFixfr1KlTzjav1xvBznq3oqIi+f1+jRkzRhkZGXrppZd07NgxPfroo5FuzRhz587Vq6++qjfffFOxsbHO9U+WZal///4R7s4csbGxncJMTEyMBg0aRMj5fwg6wDVsy5Yt+uijj/TRRx91CpBMxl69mTNn6vTp0/rlL3+phoYGpaamavPmzRoyZEikWzPGhVv3J0yYEDa+Zs0azZo169tvCNctTl0BAABjcUUjAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADAWQQcAABiLoAMAAIxF0AEAAMb6v1HPEG8+3UR0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = plt.hist(samples, bins=40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "386d6b13",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
